FLUID SOLVER INDEPENDENT HYBRID METHODS FOR MULTISCALE 
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Abstract . In some recent works 1111 1121 we developed a general framework for the construction of hybrid 
algorithms which are able to face efficiently the multiscale nature of some hyperbolic and kinetic problems. Here, 
at variance with respect to the previous methods, we construct a method form-fitting to any type of finite volume 
or finite difference scheme for the reduced equilibrium system. Thanks to the coupling of Monte Carlo techniques 
for the solution of the kinetic equations with macroscopic methods for the limiting fluid equations, we show how it 
is possible to solve multiscale fluid dynamic phenomena faster with respect to traditional deterministic/stochastic 
methods for the full kinetic equations. In addition, due to the hybrid nature of the schemes, the numerical solution is 
affected by less fluctuations when compared to standard Monte Carlo schemes. Applications to the Boltzmann-BGK 
equation are presented to show the performance of the new methods in comparison with classical approaches used 
in the simulation of kinetic equations. 
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1. Introduction. The classical fluid dynamic models like the Navier-Stokes or the Euler 
equations are not always satisfactory when dealing with large temperature or very low densities, 
and a more detailed analysis becomes often necessary to obtain correct values of the macroscopic 
quantities. In such cases a kinetic approach based on the Boltzmann equation [5] is used. The 
introduction of such a model is closely linked with the introduction of strong difficulties from the 
numerical and computational point of view. In fact, the system of equations to solve becomes very 
large, especially in multidimensional situations, and even with computers of the last generation 
the computational cost of a direct discretization is often prohibitive. Moreover the Boltzmann 
collision term that characterizes the kinetic equation, is very hard to treat in practice due to its 
nonlinear nature and physical properties. To these aims, probabilistic techniques such as Direct 
Simulation Monte Carlo (DSMC) are extensively used in real simulations for their great flexibility, 
capability of treating different collision terms and low computational cost compared to deterministic 
schemes for kinetic equations [3 [31 HOI [22]. On the other hand solutions are affected by large 
fluctuations and, in non stationary situations, the difficulty to compute averaged quantities leads 
to low accurate solutions or very expensive simulations. However, even in extremely rarefied 
regimes the fluid dynamic equations still furnish correct solution in regions of the domain where 
the gas is not subjected to sharp gradient. The direct consequence is that domain decomposition 
methods [21 [19l [1^ which consider the problem at different scales, fluid or kinetic, in different 
part of the computational domain, is a practical way to take advantage of the physics without 
loosing accuracy. We quote also the possibility to improve domain decomposition schemes through 
a moving boundary 9, 32 , in order to follow discontinuities and sharp gradients inside the domain, 
these methods are particularly important in the simulation of non stationary problems. Clearly, 
the exact identification of the non equilibrium zones remains an hard task to deal with and an 
open research argument. 

In some recent works we proposed an alternative approach to domain decomposition methods 
based on the use of different numerical methods on the whole computational domain [TJl [T2] . We 
mention here that similar hybrid approaches have been considered in [3 114[ I15[ \22\ 131] . Even if we 
develop our methods in the case of rarefied gas dynamics (RGD) the formulation we proposed per- 
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mits easily generalizations to others fields in which kinetic and hyperbolic multiscale phenomena 
are present. In order to use the hybrid approach here described, it is essential to identify a local 
equilibrium function, like the Maxwellian distribution for RGD, either analytically or numerically. 
This local equilibrium originates a model reduction from the microscale to the macroscale formu- 
lation and allows to ignore the details of microscopic interactions in terms of simplified equations 
which describe the equilibrium system. 

The schemes here presented represent an important improvement with respect to the schemes 
developed in [ll[|12j where the limiting equilibrium method was limited by construction to a kinetic 
scheme. In the present work we generalize the approach to make the method independent of the 
fluid solver used. We point out that this generalization is not trivial since in an arbitrary fluid 
solver we miss the kinetic information on the distribution function which is present in a standard 
kinetic scheme. The main advantage is that the method in the fluid limit degenerates into a 
standard fluid solver without any additional cost of a kinetic simulation. To our knowledge this is 
the first method which satisfy this property for RGD. 

Although we will focus, in the construction of the schemes, on the simplified BGK collision 
model, in principle the schemes can be extended to the full Boltzmann collision operator through 
the use of time relaxed Monte Carlo methods [211 [22l [231 [21] ■ The basic idea consists in solving 
the kinetic model and the macroscopic model in the entire domain, the first through Monte Carlo 
techniques which are robust in the fluid limit and the latter through a deterministic scheme and 
to consider as a solution a suitable hybrid merging of the two. A remarkable feature of the new 
method is the use of the hybrid moments to correct the stochastic moments in the pure Monte 
Carlo scheme. This is an important source of fluctuations reduction in the method. 

In addition we will show that it is not necessary to keep the number of sample particles fixed 
in the Monte Carlo scheme, instead it is sufficient to describe at the particles level only the fraction 
of the solution which is far from the thermodynamic equilibrium. The immediate consequence of 
the above observation is a potential reduction of the number of samples used in the Monte Carlo 
solution, and thereby, of computational time and fluctuations. These improvements are directly 
linked with the decrease of the local Knudsen number which is a measure of the rarefaction of 
the gas. The implementation of such methodology produces numerical schemes which, in general, 
are much faster than deterministic kinetic schemes and, for flow regimes close to the fluid limit, 
also to DSMC schemes. Moreover, thanks to the general formulation of the algorithm, a domain 
decomposition technique can be directly derived forcing the Knudsen number to zero (see |13j ) in 
some regions of the domain. 

Finally, let us observe that the method here developed is based on the classical operator 
splitting for the kinetic equation. This is essential if one want to match the fluid scheme with 
standard DSMC solvers for RGD, since the latter are based on such splitting. Even if there are 
well-known limitations of such splitting when dealing with Navier-Stokes asymptotics, namely the 
time step has to be related to the Knudsen number in order to describe correctly the Navier- 
Stokes level [33], here we don't aim at an under-resolved method at all the scales but simply at 
developing a method which is asymptotic preserving (AP) in the stiff limit 16]. On the other 
hand the advantages provided by our final method in terms of fluctuations and computational cost 
reduction are essentially independent of the small scales resolution but depend only on the Knudsen 
number. Note that in principle one can improve the method coupling the DSMC solver with a 
Navier-Stokes fluid description instead of a Euler one. This coupling however is not straightforward 
and we don't explore this direction here. 

The rest of the article is organized as follows. In Section 2 we introduce the BGK equations 
and his properties. In Section 3 we recall the general structure of the hybrid methods derived in 
[TTl [12] . Next in Section 4 the fluid solver independent hybrid scheme is described together with an 
acceleration technique and two possible ways in which the equilibrium fraction can be increased. 
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Section 5 is devoted to numerical results to compare performances respect traditional Monte Carlo 
and kinetic schemes. Some final considerations and future developments are discussed in the last 
Section. 

2. The Boltzmann-BGK Model. We consider the following Boltzmann-BGK kinetic model 
dtf{x,v,t)+vVJ{x,v,t) = -l—{Mf{x,v,t)-f{x,v,t)), (2.1) 

T[X, t) 

with the initial condition 



f{x,v,t = 0) = ,fo{x,v). (2.2) 

In (j2.ip the function f(x,v,t) is non negative and describes the time evolution of the distribution of 
particles with velocity -y S M'^ and position x E C K.^ , with d and D representing the dimension 
in velocity and physical space respectively. In this simplified model the Boltzmann collision term 
is substituted by a relaxation towards equilibrium. In the sequel we will work with nondimensional 
quantities, in that case t, the relaxation frequency, can be written as 

(2,3) 

where e{x, t) is the Knudsen number. Here we assume C — 1 [51[57]. Others choices of the relaxation 
time do not change the hybrid algorithm we will describe in next section. Observe anyway that 
the ratio of deterministic and stochastic component will be a function of the relaxation time, being 
linked to the ratio of the distribution function with respect to the Maxwellian equilibrium function 
as explained in details in the next Section. In the following, for simplicity, we will skip the space 
and time dependency of the Knudsen number thus e{x,t) — e. 

The local Maxwellian function, representing the local equilibrium, is defined by 

Mf{g, T)(x, V, t) = J^^^^ exp ~ ) , (2.4) 

where g, u, T are the density, mean velocity and temperature of the gas in the x-position and at 
time t 

fdv, u^- f vfdv, f \v- u\^fdv, (2.5) 

while the energy E is defined as 

E = ^ [ \v\'fdv. (2.6) 

Consider now the BGK equation (j2.ip and multiply it for 1, v, the so-called collision in- 

variant. By integrating in v the above quantities, the equations for the first three moments of the 
distribution function / are obtained. They describe respectively the conservations laws for mass, 
momentum and energy. Unfortunately, the system obtained through the above average in velocity 
space is not closed since it involves higher order moments of the distribution function. 

Note that, formally from (|2.ip as e ^ 0, the function / approaches the local Maxwellian. In 
this case it is possible to compute analytically the higher moments of / from g, u and T. Carrying 
on this computation we obtain the set of compressible Euler equations (see [4] for details) 
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Il + V, • {qu) - 

dgu „ , \ „ 

— +Vx-{Eu+pu)=0 

p^gT, E ^ ^gT + ^g\u\'^ 

where p is the thcrmodynamical pressure while (g) represents a tensor product. Higher order fluid 
model, like Navier-Stokes, can be derived similarly 

3. Hybrid Methods. The schemes derived in this paper are based on the same hybrid 
representation defined in [12j . Here we recall only the key points of the previous method, for 
details we remind to [T^] . 

For a fixed space point x we can interpret the distribution function as a probability density in 
the velocity space (the x-dependence is omitted) 

/ + CX3 
fiv,t)dv^l. (3.1) 

Next we recall the following definition of hybrid representation [T^ 

Definition 3.1. Given a probability density f{v,t), and a probability density Mf{v,t), called 
equilibrium density, we define w{v, t) e [0, 1] and f{v, t) > in the following way 

/("'')S«'<"'')^» (3,2) 
[ 1, f("J)> M,(v.t) 

and 

/>, t) = f{v, t) - w{v, t)Mf{v, t). (3.3) 
Thus f{v,t) can be represented as 

fiv, t) = />, t) + w{v, t)Mf{v, t). (3.4) 
If we take now j3{t) = min„{u'(u, t)} and /(w, t) = f{v, t) — (3{t)Mf{v, t) we have 

J{v,t)dv^l- Pit)- 
Let us define for (iit) ^ 1 the probability density 

/>,t) 



F{v,t) 



1 - (i{t) ■ 



The case /3{t) = 1 is trivial since it implies f{v, t) — Mf{v, t). Thus the probability density f{v, t), 
can be written as a convex combination of two probability densities in the form [21[ I22j 



/(«,<) = (1 - m)nv,t)+mMf{v,t). (3.5) 



Clearly the above representation is a particular case of p.4p . 

Now we consider the following general representation, including space dependence 

f{x,v,t) = f{x,v,t) +w{x,v,t)Mf{x,v,t), (3.6) 

where w{x,v,t) > is a function that characterizes the equilibrium fraction and f{x,v,t) the non 
equilibrium part of the distribution function. This representation in general can be obtained for 
the initial data of the kinetic equation using directly Definition 1. 

The starting point of the method is the classical operator splitting which consists in solving 
first a homogeneous relaxation step 

dtrix, V, t) = --{r{x, V, t) - M^-{x, V, t)) (3.7) 

and then a free transport equation 

dj^x, v,t)+v V,r (x, V, t) = 0. (3.8) 

In a single time step At the computation of the hybrid method derived in can be summa- 
rized as follows 

• Starting from a fmiction P'{x,v,t) ~ f{x,v,t) in the form p.6p solve the relaxation step 
p.7p either analytically or with a suitable numerical time integrator for stiff ODEs, like 
backward Euler. This originates the decomposition 

r{x, V, t + At)^ Xr{x, V, <) + (1 - X)Mf{x, V, t) 

= Xf{x, u, i) + (1 - A + \w{x, V, t))Mf{x, V, t), 

with < A = X{At/e) < 1 a scheme dependent constant such that A ^ as At/e oo. 
This decomposition can be cast again in the form (|3.6p taking /""(a;, v,t + At) = Xf{x, v, t) 
and w^{x, v,t + At) = 1 — A + Xwi^x, v, t). 

1. The new value w^'{x,v,t + At) follows directly from the choice of A, so from the time 
solver used for (|3.7p . 

2. The new value f^{x, v, t+At) is computed by a Monte Carlo method simply discarding 
a fraction of the samples since < A < 1 and so 'w^{x, v,t + At) > ■w{x, v, t). 

• Starting from the function /'^(x, t) — /' (a;, v,t + At) in the form ()3.6p computed above 
solve the transport step (|3.8p . 

1. Transport the particle fraction f'^{x,v,t) by simple particles shifts. 

2. Transport the deterministic fraction w'^{x,v,t)Mf{x,v,t) by a deterministic scheme. 

3. Project the computed hybrid solution f{x, v, t + At) to the form (|3.6p using Definition 
1. 

Clearly point 3 of the transport step is crucial for the details of the hybrid method. Note that 
point 2 of the transport step involves the solution of a so-called kinetic scheme for the Euler 
equations[51[^. 

In the sequel we will describe the Fluid Solver Independent (FSI) schemes which remove the 
limitations given by the use of a kinetic scheme. One major difference with respect to the hybrid 
scheme described above is that a common value for the equilibrium fraction in velocity space has 
to be chosen l3{x,t) = min„{w(a:, t)}. 

4. Fluid Solver Independent Hybrid Methods. The key feature of FSI methods is to 
take advantage from the solution of the equilibrium part of the distribution function through a 
macroscopic scheme instead of a kinetic scheme. Besides its generality, this new feature, could, in 
principle, lead to a strong reduction of the computational time with respect to any kinetic scheme 
for the fluid equation. 
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In order to describe the FSI method we introduce the projection operator V , and, in a time step 
Ai, the relaxation operator TZ/\t and the transport operator Taj- The projection operator computes 
from the kinetic variable / (or Mf) the macroscopic averages U{x,t) — {g{x,t), pu{x,t), E{x,t)), 
thus 

r{f{x,v,t))^Uix,t), V{Mf{x,v,t)) = U{x,t), (4.1) 

since the local Maxwellian Mj has the same moments of the distribution function /. The relaxation 
and transport operators solve the relaxation and transport steps. The first has the form 

nAtifix, V, t)) = Xf{x, V, t) + {l- X)Mf{x, V, t), (4.2) 

where A — exp(— At/e), whereas the second reads 

TAtifix, V, t)) = f{x - vAt, V, t). (4.3) 

Similarly we have the approximated relaxation and transport operators TZa and 7^ . For simplicity, 
since their particular structure does not play any role in the general derivation of the method, we 
assume in the sequel that Ta — T and TZa = TZ. Note that by definition Ti-AtiMj) = Mj and so 
V{TlAt{Mf)) = V{Mf). 

4.1. A simple FSI method. Let us start from an hybrid solution in the form 

fix, V, t) = {l- P{x, i))/P(x, V, t) + /?(x, t)Mf{x, V, t), (4.4) 
where fP{x,v,t) is represented by samples so that 

Nit) 

(1 - (3{x, t))fP{x, V, t)=mPY, ^(^ - PAtMv - >^jit)), (4.5) 
where Pj (t) and (t) represent the particles position and velocity, and 

= 1^ / fix,v,0)dxdv 

is the mass of a single particle, while Mf{x, v, t) is represented analytically. Note that N{t), namely 
the total number of samples is a function of time since we keep constant during the simulation. 
This is a crucial feature of the method since if we increase P{x, t) in the representation above we 
must decrease the number of samples N{t). In practice this implies that we will not be able to 
represent exactly the fraction P{x,t) but only its approximation corresponding to integer sums of 
particles. 

Since, as described in the previous section, the first relaxation step has only the consequence 
of a change of P{x, t) in (|4.4|) we derive the method starting from the transport step. 
The transport step produces the solution 

TAtifix, V, t)) = (1 - P{x - vAt, t))fP{x - vAt, V, t) + I3{x - vAt, t)Mf{x - vAt, v, t). 

From a practical viewpoint {1 — l3{x — vAt, t))fP{x — vAt, v, t) corresponds to solve a simple particle 
shift for the Monte Carlo samples. At variance the term (3{x — vAt, t)Mf (x — vAt, v, t) corresponds 
to a Maxwellian shift analogous to that usually performed in the so called kinetic or Boltzmann 
schemes for the Euler equations [SI [28]. The hybrid solution for the moments {x, t + At) is then 
recovered as 

U"ix,t + At)=V{TAtif{x,v,t))) 

= V{TAt{{l - f3{x, t))r{x, V, t))) + V{TAt{f3{x, t)Mf{x, v, t))) 
= UP{x, t + At) + U'^ix, t + At). 
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In particular {x, t+At) corresponds exactly to the approximation of the Euler solution provided 

by a kinetic/Boltzmann scheme. 

We can state the following result (see also pU] ) 

Theorem 4.1. If we denote with {x,t + At) the solution of the Euler equations {2.1^ with 
initial data U^{x,t) — V{P{x,t)Mf{x,v,t)) we have 

U^{x, t + At)^ U^{x, t + At) + 0{At^). (4.6) 

Proof. In a time step At we can write for the Euler solution 

f/^(x, t + At) = U^{x, t) + AtdtU^{x, t) + ]^{AtfdttU^(x, t) + 0{At^) 
and similarly 

[/^(x, t + At)^ U^'ix, t) + AtdtU'^ix, t) + i(At)29ttC/^(x, t) + 0(At3). 

Clearly the zero order terms in the expansions are the same since the initial data of the Euler 
equations is simply the projection of the initial data for the transport equation 

U^{x,t) =r{(3{x,t)Mf{x,v,t)) = U^{x,t). 

Now let us consider the first order terms. We have 

dtU^ = -{V^ ■ {gu),V^ ■ {gu(g)u+p),\7:, ■ {Eu + pu)f 

and 

dtU'' = V{-v ■ V,/). 

Note that this last equation is not closed since the right hand side involves third order moments 
of /. Again, however, the two terms evaluated at the initial time t coincide since the initial 
data for the transport step is the Maxwellian fraction /3{x,t)Mf{x,v,t) and so we have the usual 
Euler closure in the kinetic term. By similar arguments one can verify that the second order 
terms evaluated at the initial time are different because of the fourth order moments appearing in 
duU^{x,t) = r{v ■ V^{v ■ V^f)). This proves gH). 
□ 

By virtue of the above result we can replace the hybrid solution for the moments after the 
transport with 

U"{x, t + At)^ UP{x, t + At) + U^{x, t + At), (4.7) 

without affecting the overall first order accuracy of the splitting method. 

This hybrid values for the moments are used to compute the new Maxwellian Mj^{x, v, t + At) 
and advance the computation. To this goal we note that the next relaxation step takes the form 

7^At(TA^(/(a;, v, t))) = \T^t{f{x, v, t)) + (1 - X)Mf{x, v, t + At) 

= A(rAt((l - (3{x, t))fP{x, V, t)) + TAt{f3{x, t)Mf{x, v, t))) 
+ (1 - X)Mf{x,v,t + At) 

= (1 - I3{x, t + At))fP{x, V, t + At)+ I3{x, t + At)Mf{x, d, t + At), 

where we set 

/3(a;,t + At) = l-A, /"(a;, , t + Ai) = TAt(/(a;, <)) (4.8) 
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with A ~ er^^l^ . This shows that in order to compute the new particle fraction we need to sample 
particles from 7At(/3(a;, t)Mf{x, v, t)). In practice this can be realized in a simple way transforming 
initially the equilibrium Maxwellian part P{x,t)Mf{x,v,t) into samples and then advecting the 
whole set of samples. 

Let us denote with T^tiPix, t)M^{x, v, t j) this set of advected equihbrium samples. Computa- 
tionally this means that at each time step me must solve the full BGK model with a Monte Carlo 
scheme (I2i together with a suitable deterministic solver for the Euler equation. We can improve 
the efficiency of the above algorithm observing that we do not need to transform into samples the 
whole Maxwellian part but only a fraction A of it, where 

A > max{A(a::,t + At)}. 

X 

As discussed before, the reason for this is that we know that at the subsequent relaxation step a 
fraction (3{x, t + At) of samples will be discarded in each cell. Thus we need only 

(1 - f3{x, t + At))TAt{P{x, t)M^{x, V, t)) = X{x, t + At)TAt(/3(x, t)Af^(x, v, t)) 

advected Maxwellian particles, which is guaranteed if in any cell before advection we have at least 
A/3(x, t)My (x, w, i) particles since 

TAt{mx.t)M'P{x,v,t))^\TAt{P{x,t)M^f{x,v,t)). 

This is of paramount importance since A vanishes as e/ At — > and so the number of samples 
effectively used by the hybrid method is a decreasing function of the ration between the Knudsen 
number and the time step. 

Starting from an initial data represented by particles a simple FSI hybrid scheme for the 
solution of the BGK equation with A constant in space and time is described in the following 
algorithm 

Algorithm 1 (FSI Hybrid Scheme). 

1. Compute the initial velocity and position of the particles {v^ ,j = 1, •., N} {p^, j = 1, . . . , iV} 
hy sampling them from initial density fo{x, v). Set mP = / / fo{x, v)dxdv/N . 

2. Given a mesh Xi, i — 1, . . . ,L with grid size Ax, and an estimate of the larger sample 
velocity Vmax = 4\/27ma^, with Tmax the maximum temperature, set At^ = Ax/vmax- 

3. Compute the initial values of the moments of the distribution function in each cell Qi, 
{gu)i, Ei, i = 1, . . . ,L. 

4-. Compute the larger time step allowed by the deterministic macroscopic scheme Ato- 
5. Set At = min{AtP,AtD)■ 
6. Setn = 0,t = 0, X = e-'^^l\ A = A, /3, = 1 - A, i = 1, . . . , L. 
7. While t < tf with t f the final chosen time. 

(a) Estimate the number of Maxwellian samples we need from XP{x,t)M'^(x,v,t). 

i. In each cell set N^^ — lromid{\(3ip^ / {mP / Ax)) and sample Nf^ equilibrium 
particles from the Maxwellian with moments {pu)'^,E^. 

(b) Perform the transport step keeping track of the particles that come from the above 
sampling. 

i. Transport all particles 

p';+' =p^ + v^At, yj. (4.9) 

a. Compute the moments Uf'"^^ and the number of particles Nf in each cell using 

only the advected particles not sampled from the Maxwellian. 
Hi. Solve the Euler equations for Uf'" = (3iU" and find Uf'^^^. 



8 



iv. Compute the new hybrid moments U"^^ — JJf'""^^ + 

(c) Perform the relaxation step. 

i. In each cell set — Iround(AA'l') and discard Nf — Nf particles, 
ii. Compute the new number of particles in non equilibrium regime, in each cell 
Nf = Nf + Alf , i — 1,...,L, with Mf the transported Maxwellian particles 
(Nf'^ transported). 
Hi. Compute the effective value Af = (7V,f + Aff )/(£<^+^Ax/mP). 
iv. Set (3, = 1- Af. 

(d) Set t — t + At. n = n + 1 and compute the updated value of At. 



Remark 1. 

• In this simple version of the FSI hybrid method the value of the equilibrium fraction 
fluctuates in each cell around the constant value P — 1 — X, thus it depends on At/e. 
We will see how to remove this limitation and make the equilibrium fraction essentially 
independent of At in the optimized version of the FSI scheme. Note that fluctuations are 
due to the fact that to have mass conservation during the relaxation step we compute the 
effective value Af and set f3i — 1 — X^. 

• To avoid bias in the algorithm we used a stochastic rounding Iround(x) of a positive real 
number x defined as 



where [x] denotes the integer part of x. 
• In the fluid limit the numerical method is characterized by the particular solver adopted 
to compute the solution of the Euler equation. Thus the order of accuracy of the limiting 
scheme is completely independent from the first order splitting procedure used to solve the 
kinetic equation. This is an advantage compared to the classical approach based on kinetic 
schemes which gives limited accuracy in time. Extensions to higher order in the non fluid 
regime are not trivial since we are limited to first order accuracy in time by Theorem 14. II 

4.1.1. Matching moments. In order to have a conservative scheme it is desirable that the 
set of advected equilibrium samples satisfies 



namely the kinetic particles solution to the fluid equations in one time step should match the 
direct solution to the limiting fluid equations. Moreover, since the right hand side is not affected 
by statistical sampling error, imposing (|4.10p will decrease the variance of the samples. 

To this goal it is natural to use a moment matching approach This can be done by 
simple transformations of the sample points. Given a set of samples i^i,...,!/,/ with first two 
moments fii and 112 and a better estimate mi and m2 of the same moments we can apply the 
transformation [31 [53] 



end while 




[x] with probability [x] + 1 — x, 
[x] + 1 with probability x — [x], 



V{TAt{l3{x, t)MP{x, V, t))) = C/^(a;, t + At), 



(4.10) 




to get 



Of course this renormalization is not possible for the mass density. In fact, to keep the algorithm 
simple, we take the weight of each particle equal and constant during the simulation and this 
implies that we can have only integer multiples of such weights as mass density values in each cell. 

However thanks to the particular structure of the algorithm we can perform also a matching 
procedure for the mass using the following trick. After the transport of Maxwellian particles we 
need in each cell, in order to perform the moment matching of order zero, a number of particles 
given by 

Mf = Iround(Ap^(xi, i + At)/mP). 

This can be done if we take A large enough before transport which guarantees that we have enough 
particles in each cell. In this way the difference between the particles mass and the Euler mass is 
below the mass of one single particle. 

Next, to have exactly mass conservation we compute the effective values 

Af = (A/f + TVf )/(p^(x„t + Ai)Ax/mf + TVf), (3^{x^,t + Ai) = 1 - Af, 

used in the method. After this we renormalize the transported equilibrium samples in each space 
cell as described above so that they have the same momentum and energy of the Euler solution. 

Similarly one can apply a moment matching strategy when sampling from the Maxwellian 
during the relaxation step. In this case, as an alternative to the moment matching technique 
described above, one can use the algorithm developed by Pullin [30] . 

Note that the whole method can be seen as a Monte Carlo scheme for the BGK equation 
in which we try to reduce fluctuations substituting the moments of the transported Maxwellian, 
computed with particles, with the moments given by the solution of the compressible Euler equation 
obtained with a deterministic macroscopic scheme. Moreover, as described above, if we force the 
equilibrium particles to follow the moments given by the fluid equations we can reinterpret the 
algorithm as a fluid-dynamic guided Monte Carlo scheme. 

4.2. Optimal FSI Methods. The method just described does not take into account the pos- 
sibility to optimize the equilibrium fraction by increasing its value in time and make it independent 
on the choice of the time step. In fact at each time step the equilibrium structure is entirely lost 
and the new fraction of equilibrium is only given by the relaxation step (see l4.8|) . However, in prin- 
ciple, it is possible to recover some information from the transported local Maxwellian although 
we know it through samples rather than analytically. We recall, in fact, that wc do not get any 
microscopic information from U^{x,t) which corresponds to the solution of the Euler equation 
with a macroscopic numerical scheme. In the sequel, we will propose a method to optimize the 
equilibrium fraction f3{x, t) after the transport step. We start describing the generalization of the 
hybrid method once this optimization has been achieved. 

Thanks to Definition 1 we can define the velocity dependent optimal equilibrium fraction as 
the ratio of the transported Maxwellian at time t respect to the new local Maxwellian at time n -I- 1 

f r..(/3(.,0M/(x,.,t))<Aff(x,.,t + At), 

w {x,v,t + At) ~ < M"(x,v,t + At) 

i 1, TAtiPix,t)Mf{x,v,t)) > Mf{x,v,t + At), 

and the optimal equilibrium fraction as 

P^ix, t + At)= miniw^ix, v, t + At)}. (4.11) 

V 

This value can be considered optimal, in the sense that it is the maximum allowed value for which 
we have a decomposition like 

TAtiPix, t)Mf{x, V, t)) = Mf{x, V, t + At)+ iS^ix, t + At)Mf{x, v, t + At) (4.12) 
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with Mf(x, v,t + At) > 0. Clearly similar decompositions hold true for any fraction of equilibrium 
below the optimal one. 

Suppose, for simplicity, that P''{x,t) = at the beginning of our computation, it follows that 
the method reads in the same way from equation (|4.4[) to equation (|4.7p . Now, given an estimation 
for f3''{x, t + At) the next relaxation step reads as 

T^At{TAt{f{x,v,t))) = XTAt{f{x,v,t)) + {l-\)Mf{x,v,t + At) 

= A(rAt((l - p(x, t))fP{x, V, t)) + TAt{f3{x, t)Mf{x, V, t))) 
+ {1- \)Mf{x,v,t + At) 

= A(rAt((l - (3ix, t))fP{x, V, t)) + P'{x, t + At)Mf{x, V, t + At) 

+ Mf{x,v,t + At)) + (1 - X)Mf{x,v,t + At) 

= (1 - (3{x, t + At))fP{x, V, t + At)+ (3{x, t + At)Mf{x, v, t + At), 

with 

Pix, t + At) = l-X{l- P^x, t + At)) (4.13) 

and 

F(x vt + At) = ^At((l - Pjx, t))p{x, V, t)) + Mfjx, V, t + At) 
^ ' ' ' 1- I3''{x,t + At) ' ^ ' ' 

In order to sample from the distribution Mf{x, v,t + At), which is obtained as a difference of 
two distribution functions, see (|4.12[) . we can sample particles, exactly as in the previous section, 
from the transported Maxwellian and then apply an acceptance rejection technique that reads 

Algorithm 2 (Acceptance- Rejection Sampling) . 

do i = \, N with N number of particles to be sampled 

1. Select randomly one particle from the distribution TAt{P{x,t)AIf{x,v,t)); 

P^ix, t + At)Mf{x, V, t + At) 

2. with probability 1 — ; — keep the particle. 

TAt[p{x,t)Mf[x,v,t)) 

In the above algorithm particles can be taken more than once, in other words the sampling is 
not exclusive. Finally a moment matching strategy similar to the one described in the previous 
section can be used in such a way that the equation 

V{Mf{x, V, t + At)) = U^{x, t + At)- f]%x, t + At)U"{x, t + At), (4.15) 

is satisfied exactly. 

The major problem we have to face when practically evaluating /3'^{x, t+At) is that {x, v, t+ 
At) is known analytically while TAt{P{x,t)Mf[x,v,t)) is known only through samples. From a 
numerical viewpoint when approximating P'^{x,t + At) we want to avoid overestimates since these 
may produce unphysical solutions. In the following description, to simplify notations, we restrict 
to 1-D in velocity and physical space, extensions of the methods to multidimensional cases are 
straightforward. Our goal is to find an estimation of j3'^{x,t) given by (14. lip . Without loss of 
generality we assume that at each point x there exist a velocity v such that 

TAt{f3{x,t)Mf{x,v,t)) < Mf{x,v,t + At). 

In fact, for those space points x where the above assumption is not satisfied we simply have 
(3''{x,t + At) ^ 1. 
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The first and simplest method consists in measuring the departure from equilibrium recon- 
structing the transported Maxwellian from samples. In order to do that we need a grid in velocity 
space and a loop over the particles inside each spatial cell. We omit here the details of the differ- 
ent reconstruction methods that can be used, we refer to [26j (and the references therein) for the 
technical aspects. 

Once we have reconstructed T/^tiPixi, t)Mf{xi, v, t)), with {xi}i^i a mesh of the physical space, 
we are able to determine the quantity 

(3 {x„t + At) = mmi — >■ (4.16) 

f Mj'(Xi,v,t + At) J 

This method presents several drawbacks. The reconstruction of the distribution function from 
samples increase the computational cost, moreover a small number of particles inside a cell, which 
is quite common in applications, gives large fluctuations and this turns in an imprecise estimate 
of P^ix^t). 

A better way to estimate the equilibrium fraction /3^(x,t) after the transport is based on the 
analysis of a deterministic transport of the Maxwellian part. Again we introduce a grid in space. 
We consider the following scheme for the transport of the Maxwellian fraction 

— [-v — — = 0, w > 0, 

At Ax , - > 

(4.17) 

— +v^^ — = 0, v<0, 

At Ax 

where M^^{v) « Mf{xi,v,t"), Mj'l^{v) « M/(xi, w, t"+^) and Ax is the mesh size in space. We 
put an hat on the transported Maxwellian to distinguish it with respect to the local Maxwellian 
at time t + At, which is accordingly to the notations, M^f^{v). The scheme described above is 
a simple first order upwind for the Maxwellian transport. Of course to effectively perform the 
computation it is necessary to truncate the Maxwellian in order to obtain finite values for the 
velocity and time step larger than zero. Typically this truncation leads to several problems which 
are common in numerical methods for kinetic equations (see for example [181 125] ) . Here we are 
only interested to estimate the departure from the equilibrium of the transported Maxwellian and 
for that scope we choose a bound for the velocity space in such a way that no additional time step 
restrictions are imposed. Solving Eqs. (|4.18p we obtain 



' . X / vAt\ „ , , vAt „ , , 



>0, 



M-f(v) = (l + ^) MUv) - ^M,V,(z;), . < 0. 



(4.18) 



Note that, since \v I At < Ax, the updated fmiction M].'+\u) is a convex combination of the local 
Maxwellian in the cells i and i — \ for positive velocities and in the cells i and i + 1 for negative 
velocities. 

Now, ignoring the error introduced by the truncation in velocity, in each cell the equilibrium 
fraction /3^(xi, t 4- At) satisfies 
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In the general case a numerical method is required to compute the minimum on the right hand side. 
This operation can be expensive since it has to be done at each time step and in each spatial cell. 
However we can restrict ourselves to a lower estimate of (3'^{xi,t + At) (to avoid an overestimate 
of the equilibrium fraction) and we can choose instead of the minimum a lower bound for (|4.19p 
using the convexity property of the scheme. That value can be estimated observing that 



1 > min 1 


^ ■ i 

, mm < 




ti>0 1 







mm < — 77 — — > > mm < mm < — 77 — — > , mm 



= /3^(x,t + At), (4.20) 



mm < — 77 — —. > > mm < mm < — „ , .. > , mm < — — —. > > 

-<o \ Mff+\v) J - \v<o\ Mff+^ (v) J ' -<o [ Mff+^ (v) J J 

= Plix,t + At), (4.21) 

and setting 

(i^ix, t + At)^ min{/?^(a;, t + At), /3£(x, t + At)} (4.22) 

where the minimum of the ratios in (|4.20p - (|4.2ip can be computed exactly, being the ratios of 
Maxwellian functions. 

An algorithm that can be used to implement the optimized FSI method for the solution of the 
BGK equation, in which for simplicity A is constant, is the following 
Algorithm 3 (Optimized FSI Hybrid Scheme). 

1. Compute the initial velocity and position of the particles = 1, •■, N} {p^, j — 1, . . . , N} 
by sampling them from initial density Jq{x,v). Set vnP ~ J J fQ{x,v)dxdv/N . 

2. Given a mesh Xi, i — 1, . . . ,L with grid size Ax, and an estimate of the larger sample 
velocity Vmax — 4^2Tma£C7 with T^ax the maximum temperature, set At^ = Ax/vmax- 

3. Compute the initial values of the moments of the distribution function in each cell Qi, 
{Qu)i, Ei, i = 1, . . . ,L. 

4-. Compute the larger time step allowed by the deterministic macroscopic scheme At^. 
5. Set At ^ min{AtP,AtD)■ 
6. Setn = 0,t^O, X = e~^*/^, A = A, A = 1 - A, i = 1, . . . , L. 
7. While t < tf with t f the final chosen time. 

(a) Estimate the number of Maxwellian samples we need from XP{x,t)Mj{x,v,t). 

i. In each cell set N^^ — h:o^md{(3ip^ / {mP / Ax)) and sample N^'^ equilibrium par- 
ticles from the Maxwellian with moments {pu)^,E". 

(b) Perform the transport step keeping track of the particles that come from the above 
sampling. 

i. Transport all particles 

p]+' =p] + i^fAt, Vj. (4.23) 

ii. Compute the moments Uf'"^^ and the number of particles Nf in each cell using 

only the advected particles not sampled from the Maxwellian. 
Hi. Solve the Euler equations for Uf'" = (3iU" and find Uf'""'^^ . 
iv. Compute the new hybrid moments U^i^^ — Uf'"^^ + jjE,7i+i 

(c) Compute the optimal equilibrium fraction [3^ '"^^ as described in ^.20\ )- [4-2^ . 
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(d) Perform the relaxation step. 

i. In each cell set = Iround(AA^f ) and discard Nf — particles, 
ii. In each cell sample N^ '^ particles from the distribution AI'^^^{v)with the acceptance- 
rejection technique described in Algorithm 2. 
Hi. Apply the moment matching technique to the N^'^ particles in order to satisfy 

iv. Compute the new number of particles in non equilibrium regime, in each cell 

Nf = Nf + N^^,i^l,...,L. 
V. Compute the effective equilibrium fraction (3i — 1 — {N^ + N^'^)/{g"^^Ax/mP). 

(e) Set t = t + At, n = n + 1 and compute the updated value of At. 
end while 

Remark 2. 

• We emphasize that the first order upwind method to compute the deterministic transport 
of MaxweUians is never used in practice. It serves us only as an approximation strategy 
in order to compute a lower bound for the optimal equilibrium fraction (j'^{x,t + At). In 
this sense it is worth to notice that the additional first order dissipation introduced by 
this upwinding produces additional smearing and in principle, close to discontinuities, can 
produce overestimates of the equilibrium fraction when computed from (|4.19p . Besides 
computational efficiency this is an additional motivation to use a lower bound for that 
value. 

• The hybrid composition of the solution in the final method does not depend on the time 
step At but only on e. Note, however, that small time steps, below the CFL condition 
of the deterministic Euler solver, may increase the computational cost. To reduce this 
effect one can use different time steps in the kinetic and the Euler solver and perform the 
hybridization and matching only at intermediate steps. This strategy can be used, for 
example, where there is the need to resolve small scales at the Navier-Stokes level or in 
boundary layer effects. 

5. Implementation and numerical tests. In principle any finite volume or finite difference 
numerical scheme can be used to solve the compressible Euler equations in our hybrid method. 
In the sequel we will use a second order finite volume MUSCL type relaxed scheme (see [17] for 
details). 

In the next paragraphs we analyze the performances of the fiuid solver independent hybrid 
schemes in comparison with a classical Monte Carlo method (MCM) for several one-dimensional 
problems with different Knudsen numbers ranging from e = 10~^ to £ = 10~^. 

As a reference solution we used a deterministic discrete velocity model (DVM) for the BGK 
equation for all tests (see [18] for details). We use the shorthand FSI, FSIl and TVD to denote the 
simple FSI method, the optimal FSI method and the second order in space MUSCL Euler solver 
respectively. 

5.1. Accuracy test. Because our aim is to compare the differences in accuracy and com- 
putational time between the different methods first we have considered an accuracy test with a 
a periodic smooth solution. We compare the results of the kinetic-solver based hybrid methods 
developed in [12] to the new independent fluid solver schemes. 

We report the total Li norm of the errors for the conserved quantity p, u, and T as the 
computational times by considering a problem with the following initial data 

'2i7TX 

g{x, 0) = 1 + ag sin — — 

M(a;, 0) = 1.5 + a-u sin — — (5.1) 
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77^ — V — 


e = lO-"* 


e = 10~* 


£ = 10^^ 


MCM N=1500 


23 sec 


25 sec 


27 sec 


26 sec 


BHM N=1500 


35 sec 


25 sec 


22 sec 


22 sec 


BHMl N=1500 


34 sec 


20 sec 


19 sec 


20 sec 


BCHM N=1500 


15 sec 


11 sec 


17 sec 


21 sec 


FSI N=1500 


25 sec 


22 sec 


3 sec 


0.6 sec 


FSIl N=1500 


18 sec 


17 sec 


2 sec 


0.6 sec 


FSI N=500 


9 sec 


8 sec 


0.4 sec 


0.3 sec 


FSIl N=500 


7 sec 


6 sec 


0.4 sec 


0.3 sec 



Table 5.1 

Accuracy test. Computational times for FSI and FSIl with two different initial numbers of particles N = 1500 
and JVi = 500 compared to the previous hybrid methods (see flSt/ ) for different values of the Knudsen number. 





£ = IQ-'^ 


£ = 10-^ 


e = 5 X 10-" 


e= 10-4 


MCM 


5.494 


5.786 


5.153 


5.184 


FSI 


5.545 


3.926 


3.067 


0.268 


FSIl 


4.588 


3.406 


2.451 


0.243 



Table 5.2 

Accuracy test. Li norm of the error (in units of 10-^J for the density with respect to different values of the 
Knudsen number e. 



E{x, 0) = 2.5 + QT sin — — 

where we set 

Qg = 0.3 a„ ~ 0.1 as — 1. 

The equations are integrated for t £ [0, 5 x 10~^] using 200 space cell. 

In order to make a fair comparison with the previous schemes named BHM, BHMl, BCHM 
(see [12] for details and parameters settings of the methods) we use at the beginning 1500 particles 
for cell with the same time step of the Boltzmann-BGK schemes. Then, because in general FSI type 
schemes allow larger time steps and the moment matching techniques produce lower fluctuations 
we repeat the computation with N = 500 and the time step prescribed by the Monte Carlo method. 
We remark that, while estimating the time differences between the FSI schemes and MCM is quite 
easy (in fact the FSI methods are based on the MCM for the kinetic part), the same comparison 
with another kinetic solver such as DVM or BHM, is not straightforward due to the several possible 
choices involved in such schemes (for example the way the velocity domain is truncated and the 
type of solver chosen for the space derivatives). For this reason we stress that the simulations 
times, reported in Table 1, are just indicative. It is clear that the hybrid methods here developed 
represent a strong improvement with respect to the previous schemes as well as to the classical 
Monte Carlo method in terms of computational time. In Figure [57T] we report the total number of 
particles used by the different algorithms. Note how computational time and fluctuations reduce 
dramatically when the Knudsen number diminishes. 

The results for the relative Li errors are reported in Tables !??^ 15. 31 respectively for density, 
mean velocity and temperature for the FSI, FSIl and the Monte Carlo schemes. In all the methods 
we used N = 200 particles for cell and the moment matching techniques. We notice that the hybrid 
methods have approximately the same accuracy of the Euler solver for small Knudsen numbers 
and the same accuracy of the Monte Carlo method for large Knudsen while to intermediate values 
correspond intermediate behaviors. From the above results it is clear how the performances of 
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Number of particles in the domain lor Knudsen=0.D1 




Number of particles in the domain for Knudsen=0.0D1 




D 0.005 D.Ol D.015 D.02 0.025 0.03 0.035 0.04 0.045 0.05 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 



Number ol particles in the domain (or Knudsen=0.00D5 



Number ol particles in the domain for Knudsen=0.00D1 



i«i xoxoccccoooooxoooaxoxico(x^^ 




0.005 D.Ol 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 



Fig. 5.1. Accuracy test. Number of particles in time inside the computational domain for FSIl, FSI and MC 
schemes. Knudsen numbers e = 10"""^ (top left) e = 10~^ (top right) e = 5 x 10~^ (bottom left) and e = 10~'^ 
(bottom right). 





e = 10-'-' 


e = 10-^ 


e = 5 X 10"* 


e = 10-* 


MCM 


6.565 


5.437 


5.338 


6.035 


FSI 


4.802 


4.401 


3.264 


0.641 


FSIl 


5.135 


4.102 


2.848 


0.610 



Table 5.3 

Accuracy test. Li norm of the error (in units of 10~^) for the mean velocity with respect to different values 
of the Knudsen number e. 



the hybrid schemes are better than the ones of a Monte Carlo method, moreover FSIl gives in 
general better results respect to FSI both in term of computational time and accuracy for almost 
all regimes. 

5.2. 1-D Unsteady shock. Next we consider an unsteady shock that propagates from left to 
right. The shock is produced miming a specular wall on the left boundary, thus for the stochastic 
component at each time step, particles which escape from the computational domain on the left 
side are put back in the first cell with opposite velocity and opposite position respect to zero. On 
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e = 10-^ 


£ = 10-^ 


e = 5 X 10-* 


e= 10-* 


MCM 


6.762 


7.611 


7.578 


7.316 


FSI 


7.007 


6.022 


4.500 


0.641 


FSIl 


6.662 


4.939 


3.773 


0.598 



Table 5.4 

Accuracy test. Li norm of the error (in units of 10~^j for the temperature with respect to different values of 
the Knudsen number e. 



the other side particles are injected with the initial mean velocity and temperature in a number 
which corresponds to the initial density. For the macroscopic part the usual specular and inflow 
boundary condition are used. At the beginning the flow is uniform with mass £» = 1, mean velocity 
M = — 1 and energy E = 2.5. The computations are stopped when t = 0.065, the number of cells 
are 200 in space, while the initial number of particle are 500 for each space cell. In each Figure the 
solution computed with the Euler scheme and the one computed with the DVM scheme is reported. 
The FSI, FSIl and MCM are respectively depicted for density, mean velocity and temperature. 
We observe that for large Knudsen numbers FSI (Figure 15.21 left) and FSIl (Figure 15.21 right) 
provide a small improvement with respect to MCM (Figure [5TB1 left) in term of fluctuations. When 
£ decreases the non-equflibrium part becomes smaller and both FSI and FSIl (Figure [Ol and [5^ 
contain less fluctuations than MCM. Note that, since the time step here is 0(£), the reduction of 
fluctuations in FSI scheme is essentially the same for e = 0.001 and e = 0.0005 whereas for FSIl 
scheme the solution shows a remarkable improvement as e diminishes. Finally for e = 10^* (Figure 
15. 5p we are in an under-resolved regime and both hybrid methods yield similar solutions at the 
same computational time. 

5.3. 1-D Lax Shock Tube Test. Finally we consider a Lax shock tube test with initial 
values 

/ 0.445 \ / \ 

Ul = 0.598 ,if0<a;<0.5 nR={ ,if0.5<x<l. 
\ 3.5 / \ 0.48 / 

The solution is computed with 200 grid points in space, the final time is t = 0.05. The initial 
number of particle is 500 for each space cell. Each Figure contains the DVM solution and the 
Euler solution as reference. Similar considerations to those of the previous section hold for this 
test case. Thus for large Knudsen numbers the solutions computed with the hybrid methods 
(Figure [5?8ll5.9p show small improvements compared to the Monte Carlo scheme (Figure [5.12|) . On 
the other hand when the Knudsen number becomes smaller FSI and FSIl schemes (Figure 15.101 
15. lip give a considerable reduction of fluctuations. This is especially true for FSIl method which 
demonstrates the importance of a good estimate of the equilibrium fraction /3'^ after the transport. 

6. Conclusion. In this paper we have extended the hybrid kinetic methods developed in [12] 
to the case of an arbitrary fluid solver for the equilibrium component. Although, the simplified 
BGK collision operator has been used to develop the schemes, extensions to the full Boltzmann 
operator of rarefied gas dynamics in principle are possible through the use of time relaxed methods 
pn I24j . We plan to extend the schemes to the full Boltzmann equation in the nearby future. 

The results obtained are very promising in terms of computational cost with respect to tradi- 
tional deterministic methods for kinetic equations like discrete velocity model or spectral schemes 
[TSt [25] . In addition the FSI hybrid algorithms yield less fiuctuations with respect to direct simu- 
lation Monte Carlo methods and, close to the fluid regime, they permit to compute results faster. 
A remarkable feature of the FSIl scheme is that the equilibrium fraction is essentially independent 
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of the choice of the time step and thus provides more accurate results then Monte Carlo methods 
even in resolved regimes. 

Some open questions remain on alternative ways to estimate and increase the fraction of 
equilibrium in each space cell without increasing the computational cost. It is also interesting, and 
will be the subject of future works, to measure the response of the FSI hybrid methods in others 
situations such as simulations of nanosystem devices, plasma physics problems or turbulence. 

Acknowledgements. The authors would like to thank Russ Caflisch and Pierre Degond for 
several stimulating discussions. 

REFERENCES 

[1] G. A. Bird, Molecular gas dynamics and direct simulation of gas flows, Clarendon Press, Oxford 
(1994). 

[2] J. F. BoURGAT, P. LeTallec, B. Perthame, and Y. Qiu, Coupling Boltzmann and Euler equations 

without overlapping, in Domain Decomposition Methods in Science and Engineering, Contemp. 

Matli. 157, AMS, Providence, RI, (1994), pp. 377-398. 
[3] R. E. Caflisch, Monte Carlo and Quasi-Monte Carlo Methods, Acta Numcrica (1998) pp. 1-49. 
[4] C. Cercignani, The Boltzmann Equation and Its Applications, Springer- Vcrlag, New York, (1988). 
[5] C. Cercignani, Rarefied Gas Dynamics: From Basic Concepts to Actual Calculations, Cambridge 

Texts in Applied Mathematics, Cambridge University Press, Cambridge (2000). 
[6] F. CORON, B. Perthame, Numerical passage from kinetic to fluid equations, SIAM J. Numer. Anal., 

vol. 28 (1991), pp. 26-42. 

[7] N. Crouseilles, P. Degond, M. Lemou, A hybrid kinetic-fluid model for solving the gas-dynamics 
Boltzmann BGK equation, Journal of Computational Physics, vol. 199 (2004), pp. 776-808. 

[8] S. Deshpande, a second order accurate kinetic theory based method for inviscid compressible flow. 

Journal of Computational Physics, (1979). 
[9] P. Degond, G. Dimarco, L. Mieussens, A moving interface method for dynamic kinetic-fluid 
coupling. Journal of Computational Physics Vol. 227, pp. 1176-1208. 

[10] P. Degond, S. Jin, L. Mieussens, A Smooth Transition Between Kinetic and Hydrodynamic Equa- 
tions , Journal of Computational Physics, vol. 209 (2005), pp. 665-694. 

[11] G. Dimarco, L. Paresciii, Hybrid multiscale methods I. Hyperbolic Relaxation Problems, Comm. 
Math. Sci., 1, (2006), pp. 155-177. 

[12] G. Dimarco, L. Pareschi, Hybrid multiscale methods 77. Kinetic equations, SIAM Multiscale Mod- 
cling and Simulation Vol 6., No 4,pp. 1169-1197, (2008). 

[13] G. Dimarco, L. Pareschi, Domain decomposition techniques and hybrid multiscale methods for ki- 
netic equations, Proceedings of the 11th International Conference on Hyperbolic problems: The- 
ory, Numerics, Applications, pp. 457-464. 

[14] W. E, B. Engquist, The heterogeneous multiscale methods. Comm. Math. Sci., vol. 1 (2003), pp. 87- 
133. 

[15] W. E, B. Engquist, Multiscale Modeling and Computation, Notices of the AMS, vol. 50(9) (2003), 
pp. 1062-1070. 

[16] S. Jin, Runge-Kutta methods for hyperbolic conservation laws with stiff relaxation terms, J. Comput. 

Phys., 122 (1995), pp. 51-67. 
[17] S. Jin, Z. P. Xin, Relaxation schemes for systems of conservation laws in arbitrary space dimensions. 

Comm. Pure Appl. Math., vol. 48 (1995), pp. 235-276. 
[18] L. Mieussens, Discrete Velocity Model and Implicit Scheme for the BGK Equation of Rarefied Gas 

Dynamic ,Mathematical Models and Methods in Applied Sciences, Vol. 10 No. 8 (2000), 1121-1149 
[19] P. LeTallec, F. Mallinger, Coupling Boltzmann and Navier-Stokes by half fluxes Journal of 

Computational Physics, vol .136 (1997), pp. 51-67. 
[20] S. Liu, Monte Carlo strategies in scientific computing. Springer, (2004). 

[21] L. Pareschi, R. E. Caflisch, Implicit Monte Carlo methods for rarefied gas dynamics I: The space 
homogeneous case, J. Comput. Phys., vol. 154 (1999), pp. 90-116. 

18 



L. Pareschi, R. E. Caflisch, Towards an hybrid method for rarefied gas dynamics, IMA Vol. App. 

Math., vol. 135 (2004), pp. 57-73. 
L. Pareschi, S. Trazzi, Numerical Solution of the Boltzmann equation by time relaxed Monte Carlo 

(TRMC) methods, International Journal for Numerical Methods in Fluids, vol. 48 (2005), pp. 

947-983. 

L. Pareschi, G. Russo, Time Relaxed Monte Carlo methods for the Boltzmann equation, SIAM J. 

Sci. Comput. 23 (2001), pp. 1253-1273. 
L. Pareschi, G. Russo, Numerical solution of the Boltzmann equation. I. Spectrally accurate ap- 
proximation of the collision operator. SIAM J. Numer. Anal. 37 (2000), no. 4, 1217-1245. 
L. Pareschi, Hybrid multiscale methods for hyperbolic and kinetic problems, Esaim Proceedings, Vol. 

15, T. Goudon, E. Sonnendrucker & D. Talay Editors, pp. 87-120, (2005). 
S. PlERACClNl, G. Plippo, Implicit- explicit schemes for BGK kinetic equations, Journal of Scientific 

Computing, Volume 32, Number 1, July 2007 , pp. 1-28. 
B. Perthame, Boltzmann type schemes for gas dynamics and the entropy property, SIAM J. Num. 

Anal, Vol. 27, (1990), pp. 1405. 
B. Perthame, Second-Order Boltzmann Schemes for Compressible Euler Equations in One and Two 

Space Dimensions, SIAM J. Num. Anal., Vol. 29, No. 1, (1992), pp. 1-19. 
D. I. PuLLiN, Generation of normal variates with given sample, J. Statist. Comput. Simulation, 9 

(1979), pp. 303-309. 

R. RovEDA, D.B. Goldstein, P.L. Varghese, Hybrid Euler/Direct Simulation Monte Carlo Cal- 
culation of Unsteady Slit Flow, AIAA J. Spacecraft Rockets, vol. 37 (2000), pp. 753-760. 

S. TiwARi, Coupling of the Boltzmann and Euler equations with automatic domain decomposition, J. 
Comput. Phys., vol. 144, 1998, 710-726. 

K. Xu, A gas-kinetic BGK scheme for the Navier-Stokes equations and its connection with artificial 
dissipation and Godunov method. J. Comput. Phys. 171 (2001), no. 1, 289-335. 



19 




20 



Density, Knudsen number =0.001 



Density. Knudsen nLmber=0.001 




21 




22 



Densily. Knudsen number=0.0001 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 



Density, Knudsen number =0.0001 



- FSI1 
DVM 

- TVD 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 



Velocity. Knudsen nLmber=0.0001 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 



Velocity, Knudsen number =0.0001 



- FSI1 
DVM 

- TVD 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 



Temperature, Knudsen number=0.0001 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 



Temperature, Knudsen number=0.0001 



- FSil 
DVM 

- TVD 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 



Fig. 5.5. Unsteady Shock: e = 10~^. Solution at t = 0.065 for FSI (left) FSIl (right). From top to bottom 
density, mean velocity and temperature. 
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